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Abstract 

The semi-confined Ledo-Paniselian (Eocene) Aquifer in Flanders, Belgium is recharged in the areas with higher topo¬ 
graphy, where it is covered by the Bartonian Clay. Initial conditions in these sediments were marine. Hydrogeochemical 
modelling of the groundwater type distribution in the aquifer with a reaction transport model considers recharge of fresh 
CaHC0 3 water by downward flow through the overlying clay. In the aquifer, in an upstream direction, progressively more 
freshened water types are found. A typical NaHC0 3 type occurs upstream of the brackish-fresh water interface. The dif¬ 
ferent stages of cation exchange produce a chromatographic sequence of groundwater types, which agrees well with obser¬ 
vations. Cation exchange processes occurring mainly in the percolated clay, result in a Na + increase, and peaks of K + , 
NH^ and Mg 2+ in the aquifer, which are spatially separated as a result of groundwater flow. Calcite equilibrium, gypsum 
dissolution in the clay and sulphate reduction in the aquifer have also been included in the modelling. 

© 2006 Elsevier Ltd. All rights reserved. 


1. Introduction 

The objective of this paper is to show the control 
of cation exchange reactions on the groundwater 
chemistry of the sandy Ledo-Paniselian Aquifer in 
Flanders (Belgium). Marine conditions (represented 
by a marine interstitial solution and marine cations, 
Na + , K + and Mg 2+ , adsorbed onto clay minerals) 
prevailing in the sediments before the last regression 
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at the end of Tertiary times, are being gradually 
expelled by infiltrating fresh water. Cation exchange 
processes by fresh water intrusion in detrital marine 
sediments have been well documented by previous 
studies in other areas, such as aquifers in coastal 
plains in the United States (Foster, 1950; Back, 
1966; Chapelle and Knobel, 1983; Appelo, 1994a), 
in polder areas in The Netherlands (Appelo and 
Willemsen, 1987; Beekman, 1991; Stuyfzand, 1993) 
and in Tertiary sediments in Flanders, Belgium 
(Walraevens, 1987, 1990). 

The displacement of interstitial sea water by 
fresh water in marine sediments results in a typi¬ 
cal NaHC0 3 type water upstream of the salt- or 
brackish-fresh water interface. Furthermore a 
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chromatographic pattern can appear upstream of 
the salt- or brackish-fresh water interface when a 
sequence of Na + , K + , Mg 2+ and finally Ca 2+ type 
waters is developed. Column experiments in the lab¬ 
oratory (Appelo et al., 1990; Beekman and Appelo, 
1990; Beekman, 1991) and short term field tests 
(Valocchi et al., 1981) have shown the existence of 
chromatographic patterns when sediments with 
high exchange capacity, equilibrated with salt or 
brackish water, are flushed with fresh water. 

In aquifers the chromatographic pattern may be 
obscured by flow conditions in the aquifer (for 
example, mixing of waters of different origin) which 
change the relative concentrations of the cations 
Ca 2+ , Mg 2+ , Na + and K + . So the chromatographic 
patterns are difficult to detect in actual conditions. 
But confirmation of this theory has been described 
in detritic aquifers in the United States (Aquia 
Aquifer in Maryland, Appelo, 1994a), Dutch pol¬ 
ders (Beekman, 1991; Stuyfzand, 1993) or marine 
sediments (alternating clayey and sandy layers) in 
Flanders, Belgium (Cardenal and Walraevens, 
1994; Walraevens and Cardenal, 1994). 

A reaction transport model, according to the ter¬ 
minology of Parkhurst and Plummer (1993), has 
been used to model geochemical reactions in the 
Ledo-Paniselian Aquifer in Flanders, Belgium, 
including advection and dispersion processes. Dur¬ 
ing recharge of the aquifer, succeeding water types 
will develop and be displaced with time because of 
changes in the exchange complex and groundwater 
motion. This non-steady state with respect to the 
chemical groundwater compositions in the Ledo- 
Paniselian Aquifer renders it impossible to use mass 
balance calculations. But the reaction transport 
approach employed permits spatial and temporal 
analysis of the reactions involved in the aquifer. 

2. Geology and hydrogeology of the Ledo-Paniselian 
Aquifer 

The survey area comprises the northern part of 
the provinces of East- and West-Flanders in Bel¬ 
gium. The Eocene marine Ledo-Paniselian sedi¬ 
ments are part of a sequence of alternating clayey 
and sandy subhorizontal Tertiary deposits that are 
gently dipping toward the NNE. The geology of 
the Tertiary surface below the Quaternary sedi¬ 
ments is shown in Fig. 1. The Ledo-Paniselian 
Sands are overlain by the heavy Bartonian Clay, 
making the aquifer semi-confined. The depth to 
the top of the aquifer increases toward the NNE. 


The Ledo-Paniselian Sands are underlain by the 
Paniselian Clay, followed by the leperian (or Ypre- 
sian) Sands and the leperian Clay. The latter is a 
stiff clay with a thickness of over 100 m, and can 
be considered as the base of the groundwater reser¬ 
voir (Walraevens, 1987). The Bartonian Clay, con¬ 
fining the Ledo-Paniselian Aquifer, is composed of 
several units, of which units a 1, a2, a3 are stiff clay, 
while Asb-a, si and s2 consist of sandy clay to 
clayey sand. Toward its southern boundary 
(Fig. 1), the Bartonian Clay is not developed over 
its full extent, so that e.g. at Ursel, it almost exclu¬ 
sively consists of the lowermost units Asb-a and 
al (Fig. 4). 

In natural flow conditions, the semi-confined 
Ledo-Paniselian Aquifer is recharged in the topo¬ 
graphically higher regions, in which recharge is 
occurring by infiltration through the Bartonian 
Clay. These recharge areas are shaded on the map 
in Fig. 1. The arrows indicate horizontal groundwa¬ 
ter flow directions. They were derived from the pie¬ 
zometric map in Fig. 2a. This map has resulted from 
a 3D mathematical model of groundwater flow in 
the reservoir (Walraevens, 1988). It agrees well with 
piezometric observations, taken at a time before 
groundwater flow was substantially influenced by 
pumping (mostly around 1920). High hydraulic 
heads were found in both main recharge areas indi¬ 
cated in Fig. 1. The cross-section A-A' cuts through 
the westerly recharge area around Ursel. In the 
recharge areas, the vertical downward flow through 
the overlying clay has relatively large Darcy veloci¬ 
ties, in the order of some tens of mm/a. From the 
recharge areas on, it is converted into a subhorizon¬ 
tal flow in the aquifer. In cross-section A-A' 
(Fig. 1), the horizontal flow within the aquifer has 
a Darcy velocity around 1.6 m/a in the recharge 
area. Outside the recharge areas, a slow vertical 
upward flow exists through the Bartonian Clay, by 
which groundwater is leaving the aquifer. As a 
result, the quantity of groundwater flowing in the 
aquifer is decreasing in the flow direction, leading 
to a gradual reduction of flow velocities as the 
groundwater travels further away from the recharge 
area. In cross-section A-A' (Fig. 1), to the north of 
the recharge area, the horizontal Darcy flow veloc¬ 
ity is around 0.4 m/a. In the northern part of the 
cross-section, it is reduced to only 0.09 m/a. 

Fig. 2b shows the calculated hydraulic head dis¬ 
tribution in the aquifer in the present situation with 
groundwater exploitation (Walraevens, 1988). It is 
observed to differ strongly from the natural condi- 
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Fig. 1. Geology of the survey area and groundwater flow. Hydrogeological cross-section of Ledo-Paniselian Aquifer. Groundwater flow 
pathways are shown, (a) The Ledo-Paniselian Aquifer is recharged by a downward flow through the Bartonian Clay, (b) The groundwater 
is mostly leaving the Ledo-Paniselian Aquifer by an upward flow through the Bartonian Clay. Darcy flow velocities (indicated by arrows in 
figure): (1) 1.6 m/a; (2) 0.4 m/a; (3) 0.09 m/a or less (derived from the mathematical modelling of the groundwater flow in Walraevens, 
1988). 
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Fig. 2a. Calculated hydraulic head in the Ledo-Paniselian Aquifer in natural flow conditions. Isolines every 1.0 m (meters below sea level 
(b.s.l.)). 



Fig. 2b. Calculated present-day hydraulic head in the Ledo-Paniselian Aquifer with pumping. Isolines every 1.0 m (meters b.s.l.). 


tions, especially around the easterly recharge area, 
where the hydraulic head dome in natural condi¬ 
tions is now converted into a depression cone, with 
minima down to —27.0 m below sea level (b.s.l.) as a 
result of heavy pumping. The vertical upward flow 
in natural conditions has been turned into a down¬ 
ward flow over almost the whole survey area, thus 
substantially increasing aquifer recharge from 
above. 


3. Hydrogeochemical evolution 

The considered Tertiary sediments have all been 
deposited in a marine environment. At the end of 
the Tertiary, the sea regressed from the area and 
the present-day topography developed due to fluvial 
erosion. Groundwater movement was induced, by 
recharge of precipitation in the higher regions. Cal- 
cite dissolved in the infiltrating water, and freshen- 
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Fig. 3. Groundwater quality distribution in the Ledo-Paniselian Aquifer. Water type codes are according to the classification system of 
Stuyfzand (1991). The rectangular shaded area indicates the region of the aquifer front which observations have been taken to be 
compared with model results. 


ing of the sediments took place, consisting of both 
dilution of the marine pore solution and cation 
exchange, the latter mainly in the Bartonian Clay 
in the recharge areas (Walraevens, 1990). The ini¬ 
tially formed groundwaters were pushed to the 
north, in the direction of groundwater movement 
in the aquifer, as more precipitation was infiltrating 
in the recharge areas. The groundwater quality dis¬ 
tribution in the Ledo-Paniselian Aquifer results 
from this evolution, whereby marine conditions 
are being expelled by the infiltrating fresh water 
(Fig. 3). In the westerly recharge area around Ursel, 
the groundwater type is CaHC0 3 (according to the 
hydrochemical classification of water types after 
Stuyfzand, 1991). Downward of flow (to the 
NNE), the MgHC0 3 type appears. More towards 
the north, the NaHC0 3 type is found, and still fur¬ 
ther northward the NaCl type. In the meantime, 
salinity increases towards the north, from fresh 
(F), fresh-brackish (F b ), brackish ( B) to brackish- 
salt ( B s ) (classification after Stuyfzand, 1991). Thus, 
in an upstream direction, progressively more fresh¬ 
ened watertypes are found. The chromatographic- 
sequence of cation exchange is expressed by the sub¬ 
sequent surplus of Na + , followed by K + and finally 
Mg 2+ , resulting in the NaHC0 3 and MgHC0 3 
water types. The groundwater leaking out of the clay 
and entering the aquifer nowadays in the recharge 
area around Ursel, contains only the Ca 2+ cation 
in appreciable quantities. 


The sequence of groundwater types observed in 
the Ledo-Paniselian Aquifer is in excellent agree¬ 
ment with the pattern of natural groundwater flow 
(Fig. 2a). Although recent groundwater exploitation 
has substantially altered the flow pattern in the 
aquifer (Fig. 2b), the regional distribution of 
groundwater chemistry is still reflecting the natural 
flow system. 

4. Exchangeable cations in the Bartonian Clay and 
exchange parameters 

The situation of the adsorbed cations on the Bar¬ 
tonian Clay minerals has been studied by Walrae¬ 
vens (1987) and Walraevens and Lebbe (1988). 
Two borehole cores, one in the recharge area at 
Ursel and the other at Assenede in the discharge 
area, were used for determining the cation exchange 
capacity (CEC) and exchangeable cations, especially 
from the clayey sediments. Some sediment samples 
were also taken from the Ledo-Paniselian Aquifer. 
The determination of CEC and the extraction of 
exchangeable cations were performed by means of 
NH 4 acetate at pH 7. As the salts soluble in water 
also dissolve in NH 4 acetate, they were determined 
separately, and the result was subtracted from the 
one in NH 4 acetate, delivering the exchangeable cat¬ 
ions. Possible dissolution of calcite in NH 4 acetate 
at pH 7 was accounted for by attributing the surplus 
of the sum of obtained exchangeable Na + , K + , 
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Mg 2+ and Ca 2+ in excess of CEC, to calcite dissolu¬ 
tion, and subtracting this surplus from adsorbed 
Ca 2+ . 

CEC for the Bartonian Clay is around 27 meq/ 
100 g; for the Ledo-Paniselian Aquifer, CEC is 
5.5 meq/100 g on average. The amounts of 
exchangeable cations in the Bartonian Clay at Ursel 
and Assenede are shown in Table 1. 


The Bartonian Clay at Ursel shows a very clear 
depletion in adsorbed Na + , and to a lesser extent, 
when considering unit al, also in K + , when com¬ 
pared to the boring of Assenede, outside the 
recharge area (Fig. 4a). Adsorbed Ca 2+ is higher 
at Ursel (Fig. 4b). The adsorbed cations confirm 
that marine cations in the Bartonian Clay have been 
leached out to a large extent in the recharge area. 


Table 1 

Amounts of exchangeable cations (in meq/100 g) in the Bartonian Clay at the recharge area in Ursel (units al and Asb-a) and Assenede 
(units a3, a2, al and Asb-a) 


URSEL 


Mean 


Range 


N 


ASSENEDE 


Mean 


Range 


N 


CaX 2 

MgX 2 

NaX 

KX 


12 

13 

0.2 

1.9 


5-18 
7-19 
0 . 1 - 0.2 
0.4—3.0 


15 

15 

15 

14 


11 

10 

3.9 

1.8 


7-16 
5-17 
0.5-9.0 
0.2-3.2 


18 

18 

18 

18 


After Walraevens (1987). N: number of clay samples. 


Exchangeable Na* 


Ursel 



Oligocene 

sands 


Bartonian clay 
(with sandy 
horizons s^ & S2) 


Exchangeable K+ 

Assenede 



Oligocene 

sands 


Bartonian clay 
(with sandy 
horizons Sy & S2) 


Ledo-Paniselian 

sands 


Fig. 4a. Exchangeable cations (Na + and K + ) in sediments of cored borings at Ursel (the recharge area) and Assenede (outside the recharge 
area, where marine influence still persists). 
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Fig. 4b. Exchangeable cations (Mg 2+ and Ca 2+ ) in sediments of cored borings at Ursel (the recharge area) and Assenede (outside the 
recharge area, where marine influence still persists). 


Calcium is the dominant adsorbed cation, while 
Na + is practically absent. Outside the recharge area, 
the marine influence still persists, expressed by a 
high concentration of NaX. 

However, the analyzed samples of the Bartonian 
Clay at Ursel do not show a depletion in adsorbed 
Mg 2+ (Fig. 4b). On the contrary, adsorbed Mg 2+ 
shows an increase at Ursel, when compared to 
Assenede (Fig. 4b). This is a surprising result, as 
the Ledo-Paniselian Aquifer in this recharge area 
contains CaHC0 3 -water, indicating that the flush¬ 
ing of the overlying clay has practically been 
accomplished. 

Pore water analyses are available in the bore hole 
drilled at Assenede. Since CEC and the different 
amounts of adsorbed cations in this boring are 
known, it is possible to estimate the exchange coef¬ 
ficients in the clayey sediments. The exchange coef¬ 
ficients have been computed with respect to Na + 
using the Gaines-Thomas convention (in Appelo 


and Postma, 1993). Selecting an arbitrary high value 
for K NaX (10 20 ), the coefficients are obtained from 
the general exchange reaction 

Na + + F 1 IX, <-> NaX + i _1 I i+ 
with 

K Na/I = ([NaX][I i+ ] 1/i )/([IX 1 ] 1/i [Na + ]) 

= (/W [I’ + ] !/i ) / [Na + ] ) 

where square brackets denote activities of solutes in 
the pore water ([Na + ] and [I 1+ ]) and of the 
exchangeable cations ([NaX] and [IX]), i is the cat¬ 
ionic charge, ft is the equivalent fraction of the ad¬ 
sorbed cation I 1+ with respect to CEC, and /? Na is 
the equivalent fraction of NaX with respect to 
CEC (where /i Na = NaX (in meq/100 g)/CEC). 

Once K Na/ i is solved, the association constant 
Kixi (which represents the reaction: I 1 * + iX _ <-> 
IX;, with X - the cation exchanger) is computed 
from: 
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Table 2 


Exchange parameters for main 

cations. All values 

in lo g K„,(r'+-X * 

- DC). 




Na + 

K + 

Ca 2+ 

Mg 2+ 

NH+ 

Bartonian clay 

20.0000“ 

20.5768 b 

41.1728 b 

40.8810 b 

20.6000° 

Ledo-Paniselian aquifer 

20.0000“ 

20.7000“ 

40.8000° 

40.6000° 

20.6000° 



20.9317 d 

41.0115 d 

41.0510 d 



a Reference value. 

b Computed from data in Walraevens (1987). Mean values from 17 determinations in CEC, amounts of adsorbed Na + , K + , Ca 2+ and 
Mg 2+ and pore water analysis in samples of Bartonian clay (members a3, a2, al and Asb-a). 
c Standard values in thermodynamic database (PHREEDA) of PHREEQM (Nienhuis et al., 1993). 

d Computed from data in Walraevens (1987). Mean values from 4 sediment samples in Ledo-Paniselian (member Le). Not used. 


Kjxi — (KNax/K-Na/i) 1 

The activities of individual cations ([I 1+ ]) have been 
determined from the pore water analysis using 
PHREEQE (Parkhurst et al., 1980). An overview 
with the estimated coefficients, which have been em¬ 
ployed in modelling of the Bartonian Clay, is given 
in Table 2. The results are in the same order of mag¬ 
nitude as published data (Appelo and Postma, 
1993). 

For the Ledo-Paniselian Aquifer only a few sed¬ 
iment samples with pore water analysis were avail¬ 
able, and all of them correspond to the top of the 
sandy formation. Although exchange coefficients 
have also been calculated and result in very similar 
values to those estimated for the Bartonian Clay, 
more accurate modelling of the aquifer was per¬ 
formed with the standard coefficients in the thermo¬ 
dynamic database (PHREEDA) of PHREEQM 
(Table 2). Indeed, by using these standard coeffi¬ 
cients, the calculated chromatographic peaks of cat¬ 
ions could better be matched with observations in 
the aquifer. Probably, the exchange coefficients 
derived from data just for the top of the aquifer 
are not representative for the whole thickness. 

5. Geochemical modelling of freshening of the 
Bartonian Clay and Ledo-Paniselian Aquifer - 
Development of a chromatographic pattern 

A geochemical mixing cell model (PHREEQM; 
Nienhuis et al., 1993) has been used in order to sim¬ 
ulate the freshening of both clayey and sandy layers. 
A flow line was selected, starting from the top of the 
Bartonian Clay in the western recharge area at 
Ursel, and directed towards the NNE within the 
aquifer. Model results were compared with observa¬ 
tions, taken from the rectangular shaded area indi¬ 
cated in Fig. 3. 


5.1. Hydraulic schematization 

The modelling was divided into two parts, because 
of differences between groundwater flow conditions 
in the Bartonian Clay compared to the Ledo-Panis¬ 
elian Aquifer. Groundwater flow in the Bartonian 
Clay is vertical downwards; in the Ledo-Paniselian 
Aquifer the flow is subhorizontal. Groundwater 
flow velocity in the clay is of the order of some tens 
of meters, while the flow path in the aquifer is tens 
of kilometres long. Outside the recharge areas, the 
flow velocity within the aquifer gradually decreases 
in the flow direction, as part of the groundwater 
leaves the aquifer by vertical upward flow through 
the clay. 

Freshening of the Bartonian Clay was simulated 
in a vertical column of 20 m length, divided in 4 cells 
of 5 m each. The successive watertypes leaking out 
of the Bartonian Clay were flushed through a hori¬ 
zontal column of the Ledo-Paniselian Aquifer, sim¬ 
ulating the subhorizontal flow in it. In a first 
approach, 20 km of flow in the aquifer were mod¬ 
elled by dividing the flow path in 20 cells of 1 km 
each. Differences in flow velocity within the aquifer 
were accounted for in a second stage by dividing the 
flow path into two distinct parts. After 12 km, the 
succeeding watertypes flowing out of the 12th cell 
were injected into a new column of another 13 km 
length, divided in 26 cells of 0.5 km each (Fig. 5). 
As such 25 km of flow in the aquifer were modelled. 

The vertical Darcy flow velocity through the Bar¬ 
tonian Clay in the recharge area was deduced by 
mathematical modelling of the groundwater flow 
(Walraevens, 1988) and found to be around 
0.036 m/a (ranging between 10 mm/a and 90 mm/ 
a). With a porosity of 0.3, the resulting pore flow 
velocity is 0.12 m/a. 

In the Ledo-Paniselian Aquifer, the vertical water 
motion is converted into a subhorizontal one. 
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Fig. 5. Modelling approach used to flush the Bartonian Clay and the Ledo-Paniselian Aquifer. See Table 3 (model boundaries) for more 
information. 


Within the aquifer, the flow velocity decreases from 
1.6 m/a to 0.4 m/a at about 15 km downstream. An 
average velocity (1 m/a, equivalent to 3.3 m/a of 
pore flow velocity) was selected to model the evolu¬ 
tion of groundwater composition in the aquifer in a 
first approach. This mean value was considered rep¬ 
resentative for the first part of the aquifer, for which 
an accurate modelling was performed. For the 
downflow part of the aquifer, it was necessary to 
estimate another mean value for the flow velocity 


because there the flow is very limited (less than 
0.1 m/a as derived from the flow model) and the for¬ 
mer average velocity was not representative. 

5.2. Chemical model 

An overview of the model boundaries is given in 
Table 3. For both clayey and sandy layers, the ini¬ 
tial pore water quality was sea water (analysis from 
Nordstrom et al., 1979) equilibrated with calcite and 


Table 3 

Model boundaries 

pH Ca 2+ Mg 2+ Na + K + NH+ HCC^ Cl SO 2 - log P c0 , 

Initial pore water 7.86 9.8 55.1 486 10.6 0.75 1.6 566 29.3 

Recharge water 7.32 2.1 4.17 —2.00 


V. 11 °C (present average temperature of the groundwater in the recharge area); concentration in mmol/kg H 2 0. 



Flow vel. a 
(m/year) 

Dispers 

(m) 

Diffusion 

(m 2 /s) 

Cell 

length 

(m) 

Total 

length 

(m) 

Time 

step 

(year) 

CEC 

(meq/100 g) 

Chemical reactions 

Bartonian 

clay 

0.12 

1.5 

IE—09 

5 

20 

41.7 

27 

0.45 mM/cell gypsum dissolution 
calcite equilibrium 

Ledo-P.aq. 
(0-20 km) 

3.35 

1500 


1000 

20,000 

299 

5.5 

0.45 mM/cell SO4- reduction b 
+ 0.9 mM/cell organic matter b 
calcite equilibrium 

Ledo-P.aq. 
(12-25 km) 

0.29 

1500 

IE-09 

500 

13,000 

1750 

5.5 

calcite equilibrium 


a Pore flow velocity. 
b Only in the first 4 cells. 
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0.75 mM/kg H 2 0 of NH^ production. The 
exchange complexes in the clay (around 27 meq/ 
100 g) and in the sandy aquifer (around 5.5 meq/ 
100 g) were equilibrated with this solution of con¬ 
stant composition. For the recharge water, pure 
water equilibrated with calcite and a P c o 2 of 
10“ 2O atm, which at 11 °C means 2.1 mM 
Ca(HC0 3 ) 2 , was used. 

5.3. Freshening of the Bartonian 

The Bartonian Clay column was flushed with the 
recharge water. A diffusion coefficient of 10 -9 nr/s 
and a dispersivity value of 1.5 m were used (taken 
at around 10% of flow path length). During the 
flushing, cation exchange and calcite equilibria were 
maintained. 

In the Bartonian Clay the presence of gypsum is 
observed, resulting from the oxidation of pyrite con¬ 
tained in the clay. This pyrite oxidation probably 
occurred mainly during glacial periods, when the 
water table was lowered (Van Camp and Walrae¬ 
vens, 1999) and conditions were more strongly oxi¬ 
dizing. Sulphate concentrations between 100 and 
200 mg/L are usually found in the Ledo-Paniselian 
Aquifer in the southern part of the survey area, so 
a limited source of gypsum dissolution has been 
maintained during simulation (0.45 mmol/kg H 2 0 


per cell). This fixes the SO^ concentration in the 
recharge water of the aquifer to 170 mg/L. 

Fig. 6 represents the chemical composition of the 
different solutions which are leaving the Bartonian 
Clay and the distribution of the exchangeable cat¬ 
ions in the clay during the flushing. After the clay 
has been flushed 75-100 times (or pore volumes), 
the distribution on the clay’s adsorption complex 
is similar to the one found in the boring at Ursel. 
As is shown in Fig. 6, there is a loss of NaX, accom¬ 
panied by an uptake of Mg 2+ on the clay minerals, 
although the diluting solution in the modelling does 
not contain Mg 2+ . This is due to charge effects: as 
the solution concentration decreases as a result of 
dilution, the sum of cations declines, while equilib¬ 
rium with the exchanger is maintained; this induces 
a shift in the concentration ratio of cations with 
unlike charge, leading to an increase of the adsorbed 
divalent ion at the expense of the monovalent ion 
(Appelo, 1994b). So the model predicts the increase 
in MgX 2 , which explains the enrichment in this 
adsorbed cation in Ursel with respect to the Assen- 
ede boring. At this stage almost all NaX is 
exhausted, while a substantial amount of KX still 
remains, again consistent with the situation 
observed at the Ursel boring. 

At that moment (after 75-100 pore volumes have 
been flushed) the pore water solution which is leav- 



H!| distribution on the Bartonian clay exchange complex, 
observed in boring at Ursel 
(recharge area ot Ledo-Paniselian aquifer) 

Fig. 6a. Calculated amounts of exchangeable cations (in meq/100 g). The shaded area indicates the approximate present situation of the 
exchange complex in the Bartonian Clay at Ursel. 
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Fig. 6b. Calculated chemical composition of the pore water (in mmol/kg H 2 0, alkalinity in meq/L) in the Bartonian Clay during the 
flushing (concentrations above 15 mmol/kg in the first flushes are not shown). 


ing the Bartonian Clay contains Ca 2+ and Mg 2+ as 
main cations, with an appreciable K + concentration 
and no significant Na + content. This water is very 
different from the one found at present in the Ursel 
recharge area of the Ledo-Paniselian Aquifer, in 
which Mg 2+ and K + are almost absent. The clay 
must be flushed 300-350 times (or pore volumes) 
to obtain a water quality similar to the one found 
in the recharge area of the aquifer (fresh and 
CaHC0 3 -type). This may be the result of preferen¬ 
tial flowpaths in the Bartonian Clay where the flow 
is faster. The presence of preferential pathways has 
already been suggested to account for the relatively 
high hydraulic conductivities of the Bartonian Clay 
(on the order of 10 <J m/s), which were deduced 
from the groundwater flow model calibration, and 
which did not agree with laboratory values (Walrae¬ 
vens, 1990). The flow line in the analyzed boring 
must be part of a slower pathway and the freshening 
here is still in progress (Walraevens and Cardenal, 
1994). 

From Fig. 6 it is clear how succeeding watertypes 
show a chromatographic pattern during the freshen¬ 
ing of the clay: in the first flushes, the NaCl-water- 
type appears, then sequentially, NaHC0 3 , 
NaKHC0 3 , MgHC0 3 , CaMgHC0 3 and finally, 
CaHC0 3 -type. The different watertypes obtained 
in this way are those observed at present in the 
Ledo-Paniselian Aquifer and other detritic aquifers 
(Aquia Aquifer in Maryland, USA - Appelo, 
1994a - and Dutch polders - Beekman, 1991), 
where cation exchange processes are the main con¬ 


trol in the groundwater quality during the flushing 
of brackish or salt water aquifers by fresh water 
intrusion. 

5.4. Freshening of the Ledo-Paniselian Aquifer 

The succeeding watertypes leaking out of the 
Bartonian Clay have been used to model the hori¬ 
zontal column of the Ledo-Paniselian Aquifer. In 
order to account for the different flow velocities in 
both columns, several water compositions were 
selected and proportionally flushed. The total time 
modelled in the aquifer was the same as for the Bar¬ 
tonian Clay. Fig. 7 illustrates this procedure. This 
figure represents the Mg 2+ and K + concentrations 
of the water leaving the clay column. The rectangles 
indicate the average Mg 2+ and K + concentrations 
leached out of the clay during succeeding periods 
of time steps (a time step or shift in a mixing cell 
model means that the solutions reside a time At in 
the cell, equilibrating with the cell matrix, and then, 
the water volume is poured into the next cell mixing 
with it; Appelo and Willemsen, 1987; Appelo and 
Postma, 1993). In the example shown in Fig. 7, 
between shifts 700-900 the water leaving the clay 
has a Mg 2+ content near 0.7 mM/kg FLO. These 
200 shifts take 8.34 ka because each shift in the mod¬ 
elling of the clay column is 41.7 a. In the Ledo- 
Paniselian Aquifer each shift takes 298.5 a, so 
8.34 ka are equivalent to 28 shifts. This means that 
a water composition leaving the clay during 200 time 
steps will be introduced in the aquifer during 28 
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Fig. 7. Calculated K + (solid line] and Mg 2+ (dashed line) contents of the pore water leaking out of the Bartonian Clay during the flushing. 
Rectangles represent the selected K + and Mg + contents used to flush the Ledo-Paniselian Aquifer during specified time steps. The same is 
done for the other elements. This allows the use of only 13 water compositions (Table 4) instead of the 1500 water compositions that are 
leaking out of the clay. 


shifts. The total time used to model the clay will be 
the same as for the aquifer. So the 1500 shifts 
employed to freshen the Bartonian Clay are equiva¬ 
lent to 210 shifts in the aquifer and instead of using 
1500 water compositions to flush the aquifer, only 
13 representative samples are employed (Table 4). 

With these succeeding water compositions, it has 
been necessary to flush the aquifer 11 times to 
obtain the peaks of the main cations at the right 
positions along the flow path in the aquifer 
(Fig. 8). The model dispersivity was adapted to fit 
the observed shape of the peaks, to a value of 
1500 m; this value is consistent with the fraction of 
the flow path length used for modelling the Barto¬ 
nian Clay. Because the average pore water flow 
velocity in the aquifer is more than 3 m/a, diffusion 
was neglected. During the flushing, cation exchange 


and calcite equilibria were maintained. The mod¬ 
elled distribution of the main cations Na + , K + , 
Mg 2+ and Ca 2+ is independent on most other pro¬ 
cesses that also influence the water chemistry in this 
aquifer (initial -Pam gypsum dissolution in the clay, 
sulphate reduction in the aquifer) and predomi¬ 
nantly related to the chromatographic separation. 
In fact, Cardenal and Walraevens (1994) obtained 
similar results (with respect to the main cations) 
using an initial P c o 2 of 10~ 2 5 and the standard 
exchange coefficients (in PHREEQM) during fresh¬ 
ening of the Bartonian Clay. Yet, calcite dissolution 
is important: cation exchange causes a second stage 
of calcite dissolution, as Ca 2+ is removed from the 
solution to become adsorbed, and the groundwater 
becomes undersaturated with respect to calcite 
(Walraevens, 1987). This is especially important in 


Table 4 

Selected water compositions of the pore water leaving the Bartonian Clay used to freshen the Ledo-Paniselian Aquifer 


T.S Bart. 

T.S. Led. 

Ca 2+ 

Mg 2+ 

Na + 

K + 

NH+ 

TIC 

SOj“ 

cr 

pH 

1 

1 

9.90 

55.26 

485.80 

10.59 

0.750 

1.52 

29.71 

565.80 

7.809 

29 

4 

0.01 

0.05 

12.02 

0.26 

0.072 

7.03 

1.80 

0.02 

9.821 

170 

24 

0.02 

0.05 

9.73 

0.30 

0.054 

5.91 

1.80 

0.00 

9.433 

70 

10 

0.07 

0.15 

7.87 

0.70 

0.088 

5.35 

1.80 

0.00 

8.804 

20 

3 

0.70 

1.37 

2.59 

1.55 

0.213 

5.04 

1.80 

0.00 

7.830 

40 

6 

1.18 

2.17 

0.34 

1.06 

0.181 

4.93 

1.80 

0.00 

7.619 

40 

6 

1.45 

2.34 

0.03 

0.48 

0.109 

4.88 

1.80 

0.00 

7.542 

130 

18 

1.82 

2.14 

0.00 

0.01 

0.004 

4.72 

1.80 

0.00 

7.338 

200 

28 

2.46 

1.46 

0.00 

0.10 

0.036 

4.82 

1.80 

0.00 

7.454 

200 

28 

3.10 

0.74 

0.00 

0.00 

0.000 

4.63 

1.80 

0.00 

7.254 

200 

28 

3.46 

0.33 

0.00 

0.00 

0.000 

4.59 

1.80 

0.00 

7.213 

300 

42 

3.60 

0.18 

0.00 

0.00 

0.000 

4.57 

1.80 

0.00 

7.199 

100 

14 

3.73 

0.04 

0.00 

0.00 

0.000 

4.56 

1.80 

0.00 

7.187 


T.S. Bart.: number of time steps (41.7 a/step) for the Bartonian Clay during which the mentioned average water composition is considered 
to have been leaking from the clay. T.S. Led.: corresponding time steps (299 a/step) for the Ledo-Paniselian Aquifer. 
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Ca 2+ 




Mg 2+ 



Distance (km) 


Na + Cl' 




Fig. 8. Ca 2+ , Mg 2+ , SO;/, K + , Na + and Cl~ contents in the groundwater of the Ledo-Paniselian Aquifer along the flow path. Crosses are 
field data; solid lines represent calculated concentrations for the first part of the aquifer (LP1; pore flow velocity 3.3 m/a); dashed lines are 
calculations for the last 12-25 km of the aquifer (LP2; pore flow velocity used of 0.3 m/a). 


the NaHC0 3 watertype, which is the main water- 
type giving evidence of cation exchange. As the pro¬ 
duced Ca 2+ is subsequently adsorbed, while Na + is 
desorbed, high Na + concentrations result in this 
way. Despite this second stage of calcite dissolution, 
Ca 2+ concentrations in groundwater remain low 
because of cation exchange. 

The distribution of alkalinity, expressed as mg/L 
of HCO/, has also been calculated (Fig. 9). Walrae¬ 
vens (1987) has shown the relation between alkalinity 
and the extent of cation exchange in the Ledo-Panis- 
elian Aquifer. In order to maintain equilibrium with 
calcite, which is present in the sediment, this mineral 
dissolves when Ca 2+ is adsorbed on the clays. But 
the alkalinity values modelled with only calcite equi¬ 
librium (200^100 mg/L) are clearly lower than those 


observed. The modelling of alkalinity is improved 
by considering S0 4 reduction in the first 4 cells of 
the aquifer column. The occurrence of S0 4 reduc¬ 
tion is supported by the detection of FLS in the 
recharge area of the aquifer (Walraevens, 1987). 
The S0 4 reduction causes most of the SO 4 to dis¬ 
appear a few km downstream (Fig. 8 ). It has been 
modelled by adding 3.6 mM/kg H 2 0 of organic 
matter, with redox state equal to zero, in the first 
4 km of the aquifer. This allows improving the mod¬ 
elling of the alkalinity distribution in the aquifer: in 
the first 4 km there is an increase from 300 to 
400 mg/L of HCO/; then the alkalinity becomes 
more stable (the slope decreases) and starting from 
9 km, in the NaHC0 3 -type dominated water, the 
alkalinity rises again. 
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Alkalinity / pH 



Fig. 9. Alkalinity (expressed as mg/L of HCO/) and pH distributions in the groundwater of the Ledo-Paniselian Aquifer along the flow 
path. Crosses (HCO/J and squares (pH) are field data; solid lines represent calculations for the first part of the aquifer (LP1; pore flow 
velocity of 3.3 m/a); dashed lines are calculations for the last 12-25 km of the aquifer (LP2; pore flow velocity of 0.3 m/a). 


The pH follows a similar trend in the Ledo- 
Paniselian Aquifer (Fig. 9). Typical NaHC0 3 -type 
waters have pH values around 8.5 (Foster, 1950). 
The pH has been adequately modelled although 
slightly higher values, near 9.0, have also been com¬ 
puted in the last 5 km. 

The NH 4 content also rises between 6 and 12 km 
(Fig. 10), clearly associated with the increased Mg 
and K concentrations. Also this raised NH) content 
is attributed by Walraevens (1987, 1990) to cation 
exchange. However, sea water, which is the source 



Fig. 10. NHf content in the groundwater of the Ledo-Paniselian 
Aquifer along the flow path. Crosses are field data; the solid line 
represents the modelling results for NHf /Ca 2+ exchange, con¬ 
sidering an initial production of NHf under marine synsedimen- 
tary conditions and subsequent adsorption in the sediment at the 
sea bottom, according to Murray et al. (1978). The dashed line 
represents modelled NHf contents when considering NHf 
production in the present fresh and reducing conditions in the 
recharge area of the Ledo-Paniselian Aquifer. 


for the exchanged Na + , K + and Mg 2+ , does not 
contain a lot of NH 4 . 

Sayles and Manheim (1975) state that diagenetic 
changes of pore waters in terrigeneous sediments on 
the seafloor are dominated by reactions with 
organic matter, such as the formation of NH 4 . 
Murray et al. (1978), describe the subsequent 
adsorption of NH 4 onto the clay surfaces 

(CH 2 0) i (NH 3 ) > .(H 3 P0 4 ) z + l/2xS0 2 - 
-► xC0 2 + 1 /2xS 2 + vNH 3 + xH 2 0 + H 3 P0 4 
(organic matter) 
and 

NH 3 + C0 2 + H 2 0 + Clay 
<-> HCO^ + NH 4 - Clay (after Murray et al., 1978) 

The adsorbed NH 4 can then later be delivered to 
the percolating fresh water, in exchange for Ca 2+ . 
With an initial NH 4 production of 0.75 mM/kg 
H 2 0, it is possible to explain the concentrations 
above 2 mg/L found in some parts of the aquifer. 
The modelled peak appears associated to the K + 
peak, just between the MgHC0 3 and NaHC0 3 - 
types. 

It is also possible to explain the NHj production 
by the reducing conditions existing in the first few 
kilometres of flow path in the aquifer. Ammonium 
can be produced when organic matter is being oxi¬ 
dized incompletely. The conditions should be suffi¬ 
ciently reducing to cause the appearance of NH 4 
and not NOf nor NO,, that originate in more oxi- 
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dizing conditions. In fact, to the south of the study 
area, where the Ledo-Paniselian Aquifer is uncon- 
hned, appropriate conditions explain the occurrence 
of high NH 4 contents (>1 mg/L) in some parts of it. 
Similarly it is possible to model N H 4 production 
associated with the reducing conditions causing 
SO 4 reduction in the first 4 km of the aquifer. If at 
the same time NH 4 is produced, advective transport 
and also adsorption onto the clay minerals could 
perhaps reproduce the observed peak. However, 
modelling this way, the NH 4 content rises immedi¬ 
ately in the first 4 km; as a result of advective trans¬ 
port, the concentration is stable until the 10 th km; 
then it decreases by adsorption onto the exchange 
complex. The observed peak is not accurately mod¬ 
elled this way and it is not possible to fit the field 
data (Fig. 10). 

Both processes can contribute to the NH/ distri¬ 
bution in the aquifer, but the observed rise of its 
concentration, associated to the K + peak, is better 
explained and more accurately modelled by initial 
NH 4 production in the sea bottom sediment in con¬ 
tact with pore sea water, subsequent adsorption 
onto the clay minerals and later exchange for 
Ca 2+ from the infiltrating fresh water. 

5.5. Freshening of the aquifer in the northern part 

Although the model seems valid for the first 12 to 
15 km, the field data cannot be matched for the last 
part of the aquifer, where all parameters increase in 
concentration. This disagreement is related to the 
decreasing flow velocity in the downstream section. 
In the northern part, water seeps upward through 
the Bartonian Clay, and flushing of the aquifer is 
less extensive than calculated by the model which 
considers a uniform flow velocity. 

A smaller flow velocity is needed to fit the field 
data in the northern part of the aquifer. A method¬ 
ology similar to the one employed in the transition 
from Bartonian Clay to Ledo-Paniselian Aquifer, 


CELL 12 



Fig. 11. Chemical composition (mmol/kg FLO) of the pore water 
leaving the 12th cell (12 km downward the flow) of the Ledo- 
Paniselian Aquifer column. Six representative analyses of these 
water compositions were selected in order to flush another 
column of the aquifer and model the northern part of the survey 
area. 

was utilized to adapt both parts of the aquifer with 
significantly different velocities. The 12th cell was 
selected as the boundary between both parts. The 
evolution of the water quality in this cell during 
freshening of the aquifer has been represented in 
Fig. 11. Six water compositions were considered as 
representative (Table 5) and were proportionally 
flushed through another column of 13 km with a cell 
length of 500 m. Cell 12 is thus cell 0 of the new col¬ 
umn. With a flow velocity of 0.09 m/a (data from 
Walraevens, 1987), equivalent to 0.29 m/a of pore 
flow velocity, each time step in the new column is 
1.75 ka. 

Summarizing, the complete modelling follows the 
steps (cfr. Fig. 5): 

First, transition from Bartonian Clay to Ledo- 
Paniselian Aquifer: 

(41.7a/shift B )/(298.5a/shift LP i) = 0.14 

So, 1500 shifts B * 0.14 => 210 shifts LP i 

Then transition from the first 12 km in the aqui¬ 
fer to the downstream part of it: 

(298.5a/shift LP i)/(1750a/ shift LP2 ) = 0.17 


Table 5 

Selected water compositions of the pore water leaving the 12th cell in the Ledo-Paniselian Aquifer 


T.S LP1. 

T.S. LP2 

Ca 2+ 

Mg 2+ 

Na + 

K + 

NH+ 

TIC 

SO;- 

CL 

PH 

6 

1 

9.78 

55.09 

485.40 

10.58 

0.750 

1.52 

29.27 

565.80 

7.811 

19 

3 

1.13 

8.44 

170.80 

3.87 

0.323 

5.90 

10.17 

170.60 

8.024 

25 

4 

0.02 

0.11 

16.90 

0.39 

0.061 

11.20 

1.97 

3.37 

9.190 

100 

17 

0.02 

0.06 

11.48 

0.28 

0.052 

9.02 

1.80 

0.00 

9.422 

50 

9 

0.02 

0.06 

10.19 

0.38 

0.057 

8.42 

1.80 

0.00 

9.155 

12 

2 

0.13 

0.30 

6.98 

1.26 

0.163 

7.80 

1.80 

0.00 

8.322 


These water compositions have been used to freshen the northern part of the aquifer (between 12 and 25 km). T.S. LP1: time steps for the 
first part of the aquifer (299 a/step). T.S. LP2: time steps (1.75 ka/step) for the northern part of the aquifer. 









304 


K Walraevens et al. / Applied Geochemistry 22 (2007) 289-305 


So, 210 shifts LP i * 0.17 => 36 shifts LP2 

Here shifts B refers to shifts in the Bartonian Clay, 
LP1 to the previous modelling of the aquifer and 
LP2 to the modelling of the last part of the aquifer. 
The total time modelled is 6.3 ka in the three simu¬ 
lations. After 36 time steps (LP2), the modelling can 
reproduce accurately the observations in the north¬ 
ern part of the survey area (Figs. 8 and 9). 

In this part of the aquifer, cation exchange is less 
important and mixed water types are found. In these 
conditions, very limited flow combined to mixing 
with fossil sea water that remains stationary to the 
north of the survey area, takes place and diffusion 
cannot be neglected. A gradual dilution of the initial 
pore sea water occurs. So mixed hydrochemical 
facies are developed and all elements, especially 
Na + , CD and SO 4 increase. 

6. Conclusions 

By means of a reaction transport model the fresh¬ 
ening of marine sediments has been studied. This 
has confirmed the presence of preferential pathways 
in the Tertiary Bartonian Clay, suggested by 
groundwater flow modelling. The Tertiary sandy 
Ledo-Paniselian Aquifer is recharged through per¬ 
colation of this clay. The modelling was adapted 
to the particular flow conditions in the sediments: 
a slow and vertical downward flow through a thin 
clay layer (20 m thick) is converted into a faster 
and subhorizontal flow through the aquifer. 

The modelling technique employed has con¬ 
firmed the control of cation exchange on water qual¬ 
ity in the Ledo-Paniselian Aquifer and the 
importance not only of Ca 2 + /Na + exchange, but 
also of Ca 2 + /K + , Ca 2 + /Mg 2+ and even 

Ca 2+ / NHj exchange. The groundwater quality dis¬ 
tribution in the aquifer is the result of a spatial sep¬ 
aration of the marine cations being sequentially 
desorbed by the recharging fresh water (CaHCOs- 
type water), revealing a complete chromatographic 
pattern. This chromatographic pattern has devel¬ 
oped mainly in the overlying Bartonian Clay during 
the vertical downward flow that takes place in the 
recharge area. In the Ledo-Paniselian Aquifer, 
exchange reactions (although less important than 
in the clay) and SO 4 reduction occur, and the cation 
peaks are broadened by dispersion. However, in the 
downstream northern parts of the aquifer where 
flow is limited, cation exchange processes are less 
effective and gradual dilution of stationary fossil 
seawater occurs. 
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